Prédire la circulation routière à Paris

Projet ML, M2 Statistiques et Machine Learning

1 Introduction : la circulation à Paris

La ville de Paris rend publiquement accès à de nombreux jeux de données sur la circulation routière de la ville. La découverte du site web Paris Data (“Paris Data” visité le 26.12.2021) nous a motivé à utiliser cette ressource dans notre projet de prédiction. La possibilité de relier les différents capteurs dans une structure de graphe nous a intéressé. D’une part, cela motive des études du réseau entier : peut-on construire un algorithme global, qui apprend les interactions dans le réseau entier ? D’e l’autre côté’autre part, l’examination de rues individuelles mène à des questions sur la corrélation locale de la circulation. Un grand boulevard se comporte-t-il comme ses avenues voisines ? Enfin, dans tout cela, nous nous poserons la question à quel horizon temporel nous réussirons à prédire la présence de voitures dans les rues de Paris.

Avant de nous lancer, nous présentons d’abord le jeu de données en détail, sa richesse mais aussi les difficultés qu’il apporte : Sa taille et les valeurs manquantes. Ayant passé une grande partie du projet sur cet aspect, nous présentons nos stratégies pour résumer et compléter ce data set.

Cela fait, nous formulons plusieurs problèmes de prédiction, à différents horizons temporels et en utilisons différentes parties du data set. (TODO: spécifier à la fin pour avoir une petite table de matières ici)

2 Préparer le jeu de données

2.1 Les données et leurs attribus

Sur le site Paris Data, on trouve énormément de jeux de données issus des activités de la ville, notamment sur le comptage routier (“Données Du Comptage Routier” visité le 26.12.2021). Ces données sont collectées en continu grâce à des boucles électromagnétiques implantées dans la chaussée sur plus de 3000 tronçons de voies. Ces boucles et bornes de comptage mesurent principalement deux choses: le nombre de voitures qui passent et l’occupation du tronçon de route. Disponibles de 2014 jusqu’à aujourd’hui, ces deux données sont mesurées au pas horaire, résultant en environ \[ 8 \text{ ans } \times 365 \text{ jours } \times 24 \text{ heures } \times 3000 \text{ capteurs } = 2.1\times10^8 \text{ observations dans le data set.} \]

Chacune de ces lignes contient les informations suivantes, qui serons nos variables explicatives \(X_i\). Les noms des variables que nous utilisons dans notre analyse sont écrits en gras tout au long de cet exposé.

• Un timestamp t_1h indiquant date et heure de la mesure.
• Le nombre de voitures ayant passé le point de comptage pendant une heure, noté q et renommé par nbCar.
• Le taux d’occupation, correspondant au temps de présence de véhicules sur la boucle en pourcentage d’une heure, noté k et renommé par rateCar.
• La variable etat_barre indique si la route où se trouve la borne de comptage est “barée,” etat_barre=2, ou “ouverte,” 1. Nous la renommons state. Malheureusement, cette variable a souvent la valeur “inconnu,” 0, ou “invalide,” 4. La différence entre 3 et 4 est un mystère et par conséquence nous nous servons peu de cette variable.
• L’identifiant iu_ac de chaque borne. Chaque ligne contient de informations sur les compteurs adjacents.
• Le libellé (libelle) du tronçon de rue où est placé le capteur, par exemple “Champs Elysées partie B.”

Les données historiques de circulation sont enregistrées sur le site web Paris Data sous le format .txt. Chaque fichier couvre toutes les bornes pour une semaine. Pour commencer, il faut donc télécharger les fichiers de toutes les semaines et les convertir en dataframes. Ce faisant, nous effectuons également un tri par identifiant de borne iu_ac.

2.2 Simplification et aggrégation

Au vue de la quantité énorme de 200 millions d’observations, nous sommes obligés de réduire la taille du dataset. La première démarche consiste à garder les variables intéressantes: t_1h, nbCar, rateCar et iu_ac.

Ayant réduit la dimension de chaque observation, la deuxième étape est la réduction du nombre de capteurs de comptage. Pour ce faire, nous utilisons un modèle simplifié des rues de Paris. Etant donné que les capteurs couvrent les grandes axes routiers, dont le périphérique, notre objectif était d’agréger les données selon ces grands axes. Afin d’identifier où passent le plus de voitures, nous calculons les moyennes de nbCar pour chaque libellé, c’est-à-dire les moyennes à travers le temps et pour toutes les bornes associées à chaque libellé. En prenant les 200 libellés dont le nombre de voitures est le plus grand (hors périphérique), nous obtenons le graphique 2.1 qui représente très bien les axes routiers auxquels on s’attendait.
Ensuite, nous ignorons les libellés et regardons exclusivement ce graphique. De manière arbitraire, nous en déduisons le graphe de la figure 2.2, un schéma très simplifié de la circulation à Paris. Grâce à notre démarche d’utiliser les rues les plus fréquentées, cette abstraction nous permet d’avoir un modèle représentatif de la circulation parisienne.

Figure 2.1: Représentation en rouge des bornes d’observation associés aux 200 libellés les plus fréquentés (hors périphérique)

Figure 2.2: Graphe simplifié. Les points noirs indiquent les intersections entre nos arêtes.

Chaque arête dans le graphe 2.2 regroupe plusieurs capteurs dont nous agrégeons les données : après une collecte minutieuse et laborieuse des identifiants de toutes les bornes sur chacune des arêtes, nous associons pour chaque heure la moyenne de nbCar et rateCar des bornes de l’arête dont ils font partie. Cela a trois avantages:

1. Prendre la moyenne évite d’avoir des données manquantes, issues entre autres de capteurs dysfonctionnels sur une certaine période. Ce n’est pas toujours le cas, comme on le verra dans la section suivante.
2. S’il y a des capteurs où les données fluctuent beaucoup en raison de mauvais positionnement ou des perturbations indésirées, la moyenne réduit cette variance.
3. Le problème est réduit à 69 arêtes représentées par 69 dataframes.

Cette réduction des données résulte en la construction des jeux de données que l’on va manipuler et en la concrétisation du problème que l’on souhaite résoudre. Elle permet également d’obtenir des objets que nos ordinateurs sont capables de stocker efficacement. Il reste pourtant un problème: la présence de données manquantes.

2.3 Complétion de valeurs manquantes

En examinant des figures basiques des taux d’occupation, nous observons régulièrement une image comme dans le graphique 2.3 Malgré l’agrégation le long des arêtes, il reste un nombre considérable de valeurs manquantes, NA, dans le dataset. Afin d’exécuter des tâches de prédiction, il nous faut cependant des données complètes. Une partie conséquente de notre travail a donc été de compléter les valeurs manquantes.

Figure 2.3: Taux d’occupation pour une journée ouvrière avec une valeur manquante

Figure 2.4: Visualisation de valeurs manquantes pour 30 arêtes (gris = NA)

En regardant le graphique 2.4, on considère qu’il y a deux types de “trous” dans le data set:

1. Des valeurs manquantes ponctuelles comme on a pu observer dans la figure 2.3 ou celles marquées en bleu dans le graphique ci-dessus. On peut deviner qu’il s’agit de simples erreurs dans la collection des données ou des travaux d maintenance par exemple.
2. Des trous plus grands, couvrant parfois plusieurs mois, comme ceux marqués en rouge dans le graphique 2.4. Il peut s’agir d’une absence étendue de plusieurs bornes, une clôture de la route en raison d’une construction, etc. Savoir si une route est effectivement inaccessible pendant de longues période serait une tâche de recherche en soi : comme mentionné plus tôt, la variable state ne contient presque aucune information. Nous remplissons alors les grands trous sans respecter d’éventuelles fermetures, ajoutant une autre simplification.

Le fait que les arêtes ont de différentes structures de valeurs manquantes peut aussi être observé en regardant la répartition du pourcentage de NA’s à travers les dataframes dans la figure 2.5.

Figure 2.5: Distribution des valeurs manquantes

Nous observons que les données de taux d’occupation sont généralement plus complètes : ceci est dû au fait qu’il y a plus de capteurs de ce type. Par ailleurs, nous avons vérifié qu’une valeur manquante de nbCar implique toujours un NA dans rateCar.

Pour remplir les trous, la première étape consiste à ajouter des variables explicatives supplémentaires. En fonction du type de trou, nous verrons qu’ils auront un impact important sur la complétion.

2.3.1 Variables supplémentaires : Proximité temporelle et locale

Pour l’instant, nous disposons des variables nbCar et rateCar pour chaque arête et chaque heure entre 2014 et 2020. Nous exploitons deux propriétés structurelles des données pour obtenir de l’information supplémentaire:

D’abord, il s’agit de séries temporelles, donc nous rajoutons les valeurs historiques de nbCar et rateCar décalées d’une heure, d’un jour et d’une semaine (nbCarLaggedHour, nbCarLaggedDay, et nbCarLaggedWeek) comme variables explicatives. L’espoir étant que la circulation reste identique à travers ces cycles. Avec l’objectif spécifique de compléter les données manquantes, il est aussi pertinent d’ajouter ces valeurs du futur (rateCarFuturHour etc.), étant donné que le modèle d’imputation a pour but d’interpoler au lieu d’extrapoler. On fait l’hypothèse que toutes ces variables décalées seront efficaces pour remplir les trous locaux : une heure ou une journée manquante devrait facilement être inférée en utilisant des données qui sont proches dans le temps.

Pour les valeurs manquantes au milieu des grands trous, nous n’avons pas accès aux heures et jours d’avant car elles sont aussi manquantes. Cependant, grâce à l’interconnexion de nos données à travers le graphe présenté plus haut, les mesures de circulation d’une arête peuvent bien être expliquées en fonction de celle de ces voisines. Dans un effort manuel, nous avons pour chaque arête \(A\) collectionné des voisins \(V_i\), dans le sens que les voitures dans \(A\) passent ensuite par l’un des \(V_i\) ou bien l’inverse. Nous faisons ici encore des choix arbitraires de qui est voisin de qui pour deux raisons. Premièrement, inférer statistiquement quelle arête influence quelle autre nécessiterait des données déjà complètes pour une régression par exemple. Deuxièmement, utiliser toutes les autres 68 arêtes au lieu de juste 2 à 6 voisins augmenterait trop le temps d’exécution de la complétion.
Donnons un exemple illustrant le choix de voisins ainsi que la dénomination des nouvelles variables. Parmi les voitures qui entrent dans Paris par l’arête “pont amont - pont austerlitz,” il y a une partie considérable qui continue tout droit sur l’arête “pont austerlitz - chatelet.” Par conséquent, nous ajoutons la variable rateCar_pontausterlitzTOchatelet au dataframe de “pont amont - pont austerlitz.”

Jusqu’ici, nous n’avons utilisé que des informations déjà présentes dans le jeu de données. Dans la suite, nous ajoutons des variables extérieures qui pourraient également expliquer la circulation routière.

2.3.2 Variables supplémentaires : Facteurs externes

Variables temporelles

Le traffic routier étant relié à l’activité humaine, nous avons ajouté de nombreuses variables temporelles (notamment grâce à la librairie lubridate).

• year, month, day, hour résultent d’une simple décomposition des timestamps t_1h,
• time est un simple compteur de journées, commençant à \(1\) le 1er janvier 2014,
• toy (time of year), il s’agit d’un numéro entre 0 et 1 indiquant la position de l’observation dans l’année en cours,
• weekdays le numéro du jour de la semaine et weekendsIndicator l’indicatrice si le jour est un jour du weekend,
• winterHolidaysIndicator et summerHolidaysIndicator les indicatrices des vacances d’hiver et d’été définies à partir de (“Les Archives Du Calendrier Scolaire” visité le 26.12.2021),
• bankHolidaysIndicator l’indicatrice des jours fériés définie à partir de (“Jours Fériés” visité le 26.12.2021).

Index de la situation sanitaire en rapport avec le Covid-19

Nous avons également ajouté, covidIndex, un index allant de 0 à 100 représentant les restrictions du gouvernement sur la situation sanitaire en rapport avec le Covid-19. Il est calculé à partir de nombreux indicateurs et est fourni par l’université d’Oxford (“Covid Index” visité le 26.12.2021). Cette variable permettra éventuellement une étude de l’année 2020 qui voit une circulation perturbé. Pour cela, il faudra utiliser des modèles qui s’adaptent très vite à l’influence de nouvelles variables.

Météo

A partir de données de l’Organisation Météorologique Mondiale (“Observations Météorologiques Historiques En France” visité le 26.12.2021), nous avons extrait 2 variables météorologiques : temperature la température en Kelvin et precipitation les précipitations dans les 3 dernières heures en mm. Ces données ont été mesurées à Athis-Mons en Essonne.

On dispose d’un relevé tous les 3 heures environ donc on procède à une interpolation pour compléter les données. Etant donné que les mesures sont uniformément réparties, on utilise une interpolation linéaire basique à l’aide de la fonction na.approx de la librairie zoo.

2.4 Complétion via Random Forest

Avec un dataframe enrichi à notre disposition pour chacune des arêtes du graphe de circulation, nous pouvons passer au remplissage des valeurs manquantes. Nous sommes pourtant restreint dans notre recherche d’un algorithme d’imputation car plusieurs variables explicatives contiennent également des NA. D’une part, nous voulons remplir deux variables en même temps (nbCar et rateCar). D’autre part, les voisins n’ont pas moins de valeurs manquantes, laissant de grands trous dans le dataframe.
Heureusement, il y a un modèle qui peut être entraîné en dépit de valeurs manquantes parmi les variables : les forêts aléatoires constituées d’arbres CART. Une règle de décision dans un tel arbre peut être ignorée si la valeur de la variable nécessaire n’est pas renseignée. Le package miceRanger profite de ce fait en utilisant un algorithme dit d’imputation multiple : en commençant par la variable \(V_1\), une forêt aléatoire est entraînée sur les lignes où \(V_1\) n’est pas manquant. Puis, les autres lignes sont “prédites” par la forêt. Avec le nouveau dataset moins vide, le processus est reproduit et ainsi de suite (Van Buuren 2018). Nous nous arêtons pourtant à la deuxième itération vu qu’il n’y a pas d’intérêt à compléter les données des voisins.

Comme nous répétons le processus de complétion 69 fois, le temps d’exécution est de quelques heures sur nos ordinateurs. Nous n’avons donc pas le luxe d’optimiser des hyperparamètres comme nous ferons plus tard pour la tâche de prédiction. Le choix de 100 arbres par forêt et 7 variables considérées à chaque split semble un bon compromis. (voir le fichier preprocessing/missing_data.R sur git pour le code exact).

Au délà des capacités de complétion, les forêts aléatoires permettent de calculer des scores d’importance pour chaque variable explicatives. Ces scores peuvent aider à identifier lesquelles des variables ont le plus contribué à l’estimation de la cible (le score est calculé comme la réduction de variance à chaque split). Pour une arête du graphe, jussieu - saint michel, nous illustrons les variables qui ont le plus contribué à l’imputation :

Figure 2.6: Variables avec le score d’importance le plus haut pour rateCar de l’arête ‘’jussieu - saint michel’’

Figure 2.7: Variables avec le score d’importance le plus haut pour nbCar de l’arête ‘’jussieu - saint michel’’

Sur l’arête dont traitent les deux plots, on observe que rateCar et nbCar exhibent deux comportement assez différents. Pour le taux d’occupation, ce sont presque exclusivement les valeurs temporellement décalées qui contribuent à l’imputaion. Comme il y seulement 0.8% de valeurs à compléter, ceci confirme notre hypothèse que les petits trous d’une ou plusieurs heures sont très bien complétés par les valeurs des heures juste avant ou après. Pour le nombre de voitures, où la proportion de NA était plus importante, les deux variables ayant le plus réduit les variances aux splits sont des variables nbCar issues d’arêtes voisines. Les trous importants en nbCar entre Jussieu et Saint Michel semblent être le mieux expliqué par ce qui se passe autour. Cela confirme notre choix de variables pour combler les 2 types de trous. Un phénomène similaire peut s’observer dans les autres arêtes.(TODO Guillaume: Quoi d’autre écrire sur ceci ? Y a des scores intéressants dans RF (genre error Out of Bag) que j’aurais dû retenir de l’imputation ?)

Après l’imputation arête par arête, nous pouvons enfin compléter toutes les données grâce au fait que les variables reliées aux voisins ont été imputées dans les dataframes qui correspondent à ces voisins. Il est important de noter ici que nous allons dans la suite couper le jeu de données en deux (train - test) et que cette imputation sera effectuée individuellement sur chacune des deux parties. Faire ceci deux fois séparément est nécessaire pour préserver l’indépendance et permet ensuite d’estimer les erreurs de nos modèles.
Dans la suite, nous supposons que les données sont complètes. Il faut pourtant remarquer que ces forêts aléatoires font indirectement partie des modèles de prédiction que nous allons utiliser. De plus, nous gardons à l’esprit que nous avons introduit un biais dans nos données d’apprentissage comme de test.

3 Problèmes et modèles de prévision

Avec un jeu de données complet, nous pouvons enfin utiliser des algorithmes de Machine Learning pour faire de la prévision. Vu la complexité du jeu de données, il n’est même pas évident quelles questions se poser. Nous allons alors commencer par la définition de différents problèmes et de métriques de comparaison avant d’entraîner des modèles.

3.1 Définition des problèmes

Il y a deux dimensions principales du problème que nous allons examiner dans la suite, en utilisant de différents modèles d’apprentissage:

1. La première dimension du problème est temporelle : Quel moment dans le futur voulons-nous prédire ? Autrement dit, jusqu’à quel horizon temporel est-ce que nous nous permettons d’utiliser le passé pour prédire le présent ?
   1. Le cas une heure correspond à un problème du genre “Google Maps” : si nous prenons la voiture maintenant, quelle sera la circulation dans les prochaines 60 minutes ?
   2. Etant donnée la circulation du jour précédent, la tâche devient intéressante pour les régulateurs des feux de circulation ou pour la police qui pourront ainsi éviter des bouchons.
   3. Si nous voulons prédire une semaine voire un mois en avance, le modèle devra comprendre davantage le comportement habituelle de la circulation routière. Ce sera plutôt un pari basé sur des fluctuations habituelles
2. La deuxième dimension est reliée au réseau routier : Pouvons-nous profiter de la structure du graphe pour améliorer la prédiction ?
   1. La question de base sera si on réussit à battre des benchmarks sur une seule arête.
   2. Puis nous pourrons reutiliser les “voisins” introduits pendant l’imputation.
   3. Finalement, une question est si nous réusissons à trouver un seul modèle décrivant toutes les arêtes en même temps, une sorte de prédiction multivariée.

Avant de nous lancer, nous précisons comment nous évaluons nos modèles et quels sont les benchmarks, les modèles naïfs, que nous voulons battre.

3.1.1 Evaluation des modèles

Train - test - split

On découpe notre jeu de données en 2 parties de proportion 2/3 et 1/3 : une partie apprentissage de 2014 à 2017 et une partie test de 2018 à 2019. Faire le découpage de cette manière assure qu’on est évalué sur les prédictions du futur plutôt que du passé. L’année 2020 est omise à cause du Covid, mais elle pourrait faire partie d’un autre projet sur le jeu de données en relation avec la pandémie !

Mesurer la performance

Pour évaluer les performances des modèles, nous considérerons le Root Mean Square Error (RMSE) qui représente la racine carrée du deuxième moment d’échantillonnage des différences entre les valeurs prédites et les valeurs observées. Une autre métrique souvent utilisée est le Mean Absolute Percentage Error (MAPE) qui n’est pas applicable à notre situation car un nombre important des valeurs de rateCar sont nulles et on ne peut pas diviser par elles.

Modèles naïfs à battre

Afin de pouvoir comparer nos modèles à une référence, une sorte de benchmark, on construit 3 modèles témoins. Ces modèles sont dit naïfs car extrêmement simpliste :

• naiveModel1 prévoit le taux d’occupation d’une heure \(t \in \{ 0, \dots, 23 \}\) d’un jour \(j \in \{ 1, \dots, 7\}\) d’un mois \(m \in \{ 1, \dots, 12 \}\) d’une année \(a \in \{ 2018, 2019 \}\) en moyennant les taux d’occupation de l’heure \(t\) du jour \(j\) du mois \(m\) pour toutes les années \(a \in \{ 2014, \dots, 2017 \}\).

• naiveModel2 fait de même en calculant la médiane et non la moyenne.

• naiveModel3 est simplement le taux d’occupation à l’heure précédente rateCar_LaggedHour.

Pour ces trois modèles simples, on calcule leurs erreurs RMSE pour chacune des 69 arêtes pour les moyenner après (voir le fichier naive_model.r pour le code). Ceci nous donne :

Table 3.1: Performances des modèles naïfs témoins

naiveModel1	naiveModel2	naiveModel3
5.844	5.899	3.955

Les 2 premiers modèles naïfs ont des scores similaires. Le troisième modèle naïf, qui consiste à prendre la variable rateCar de manière décalée, est meilleur. On interprète cette bonne performance par la très forte corrélation du taux de circulation d’une heure avec la suivante. Dans la suite, il sera intéressant de comparer ce modèle naïf avec des modèles plus complexes pour justifier leur utilisation.

Validation croisée

Pour la séléction de modèles et plus particulièrement pour la recherche de paramètres (d’un arbre CART par exemple), nous nous servirons beaucoup de la validation croisée. Ceci n’est pas évident pour les séries temporelles comme les données ne sont pas indépendant d’un ligne à l’autre. En plus, il vaut mieux éviter l’utilisation de données futurs pour prédire le passé. Entre ces problèmes et l’efficacité algorithmique, nous utilisons trois approches différentes:

• Validation croisée sur blocs regroupant des mois successifs: on divise les 4 années de données d’apprentissage en blocs de 3 mois ou 6 mois (correspondant à 16 et 8 “CV-folds” respectivement), ce qui englobe les différentes saisonnalités de nos données. C’est la méthode qui nous semble le plus pertinent pour les séries temporelles car les blocs sont à peu près indépendants entre eux et on respecte les saisonnalités dans le découpage. Pourtant, il n’y a pas de package efficace dans R pour faire ceci de manière efficace et le fait de devoir le programmer manuellement chaque fois nous empêche de l’utiliser souvent.

• Validation croisée progressive: on fixe une fenêtre initiale d’apprentissage de 2 ans, que l’on incrémente dans l’ordre chronologique de 2 mois à chaque itération, pour entraîner le modèle et mesurer sa performance à l’horizon 1 (c’est-à-dire prédire le taux d’occupation de l’heure suivante). Pour utiliser cette méthode, on utilise le package caret avec le paramètre timeSlice.

Figure 3.1: Validation croisée progressive (“Progressive Cross-Validation” visité le 22.02.2022)

• Validation croisée habituelle : on mélange toutes les lignes et les découpe en K blocs. Cette approche habituelle pour des données iid est moins pertinente pour les séries temporelles, on utilise particulièrement le future pour prédire le passé etc.

Si nous ignorons la première idée numériquement plus dure, il semble plus naturel d’utiliser l’approche progressive pour notre type de données. Il y a pourtant le souci de largement surévaluer les données des premières deux années dans le score de CV, vu qu’elles sont représentées dans chaque couche. En plus, les différentes couches (folds) seront alors fortement corrélées comme elles ont plus de la moitie de leur contenu en commun. Les deux idées ont alors leurs failles.

Après tout, lors de nos recherches de paramètres, nous avons observé que les deux méthodes sélectionnent des paramètres assez similaires. La plupart du temps, nous nous contentons alors de simplement utiliser la validation croisée habituelle avec 16 couches. Ce chiffre est inspiré par l’approche de blocs sur mois successifs. L’avantage de ceci est que le package caret permet une implémentation efficace. Nous notons ici que nous fixons toujours une seed du générateur aléatoire pour le mélange des lignes effectué par caret. Ainsi le lecteur pourra relancer le code lui-même.

3.2 Google Maps : Prédire la prochaine heure

Commençons par la tâche qui est censée être la plus simple : Prédire la circulation dans les prochaines 60 minutes.

3.2.1 GAM pour une arête

Au début, nous nous restreignons à une arête isolée dans le sens que nous ne prenons pas en compte la circulation sur ces voisins. Notre objectif initial était de faire une simple régression linéaire, mais il y a quelques relations pas linéaires du tout. Si on essaie par exemple, d’écrire le taux d’occupation d’une avenue en fonction de l’heure, le graphique suivant montre bien que les heures de pointe seraient mieux décrites par un polynôme que pour une droite :

Figure 3.2: Une journée typique sur l’arête

Pas conséquent nous travaillons plutôt sur des GAM que sur de simples régression. Constituer un modèle additif et en sélectionner le meilleur est une tâche qui peut être longue. C’est pourquoi nous adapterons une approche greedy pour la sélection de modèles où nous rajoutons une variable explicative à la fois. Idéalement ce sera celle qui apporte le plus de puissance prédictive supplémentaire (forward search).

Pour prédire une heure dans le futur, il est incontestable que l’heure juste avant contient le plus d’informations. La variable rateCar_LaggedHour sera donc la première à considérer. L’effet de l’occupation d’une heure à l’autre n’est pourtant pas homogène à travers les heures de la journée. En plus, comme nous voyons dans le plot supérieur gauche, il y a peut-être une relation linéaire entre les variables, qui a pourtant une grande variance, voir l’image à gauche en haut :

Les autres trois plots que nous venons de voir représentent la même comparaison rateCar_LaggedHours vs. rateCar, mais sur différentes heures. On a bien l’impression que des régressions linéaires pourront nous servir, vu que la variance qui est visiblement réduite. Comparer le nuage de points à 8 heures avec ceux de 12 et de 13 heures nous donne envie de regrouper certaines heures. Ce phénomène était déjà observable dans le graphique 3.2 : dans l’heure de pointe matinale (morning rushhour), la circulation augmente beaucoup entre 7 et 8 heures tandis qu’entre 12 et 13 heures elle reste plutôt constante.

Afin de regrouper les heures de la journées en périodes semblables, nous procédons de manière heuristique : pour chaque heure de la journée, nous effectuons une régression linéaire et comparons les coefficients (intercept et slope). Ceci donne le tableau suivant.

(#tab:linReg_heureclustering)LinReg coefficients pour chaque heure

hour	0.000	1.000	2.000	3.000	4.000	5.000	6.000	7.000	8.000	9.000	10.000	11.000	12.000	13.000	14.000	15.000	16.000	17.000	18.000	19.000	20.000	21.000	22.000	23.000
intercept	-0.148	0.038	-0.218	0.039	0.007	0.163	0.414	1.155	-0.630	-1.185	1.691	2.253	1.688	2.278	0.756	1.778	2.442	1.527	1.230	1.529	2.020	0.709	2.252	1.675
slope	1.050	0.749	0.736	0.651	0.645	0.578	0.575	0.568	2.993	2.298	0.846	0.888	0.949	0.879	0.808	0.842	1.046	1.013	0.916	0.810	0.795	0.685	0.274	0.438

Les résultats sont congruents avec le graphique 3.2, par exemple l’explosion de slope à 8 heures où commence l’heure de pointe matinale. En nous basant sur ce tableau, nous regroupons de manière heuristique les heures en différents groupes: “night” de 0 à 6 heures, “7heures” à huit heures, “8heures” et “9heures” pareil, “noon” de 10 à 15 heures, “afternoon” de 16 à 17 heures, “evening_rush” de 18 à 20 heures et “evening” de 21 à 23 heures. Ainsi, chaque ligne du dataframe reçoit le nouveau attribut hour_groups. Avec ces groupes, nous pouvons établir notre premier modèle mod1, une régression linéaire avec un “mixed effect” hour_groups * rateCar_LaggedHour :

cluster_hours = function(h){
  if (h %in% c(0,1,2,3,4,5,6,7)){return("night")}
  if (h %in% c(7)){return("7heures")}
  if (h %in% c(8)){return("8heures")}
  if (h %in% c(9)){return("9heures")}
  if (h %in% c(10,11,12,13,14,15)){return("noon")}
  if (h %in% c(16,17)){return("afternoon")}
  if (h %in% c(18,19,20)){return("evening_rush")}
  if (h %in% c(21,22,23)){return("evening")}
  else{return("asfd")}
}
df_train$hour_groups = as.factor(sapply(as.factor(df_train$hour), cluster_hours))
df_test$hour_groups = as.factor(sapply(as.factor(df_test$hour), cluster_hours))
# for (group in unique(df_train$hour_groups)){
#   indices = df_train$hour_groups == group
#   plot(df_train[indices,]$rateCar_LaggedHour, df_train[indices,]$rateCar,
#        main = group, xlab = "rateCar_laggedHour", ylab = "rateCar")
# }
form1 = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour
mod1 = lm(data = df_train, formula = form1)
summary(mod1)

## 
## Call:
## lm(formula = form1, data = df_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.480  -0.836  -0.301   0.749  51.016 
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                0.846441   0.024047   35.20   <2e-16
## hour_groups8heures:rateCar_LaggedHour      2.201884   0.038810   56.73   <2e-16
## hour_groups9heures:rateCar_LaggedHour      1.935035   0.013804  140.18   <2e-16
## hour_groupsafternoon:rateCar_LaggedHour    1.097172   0.003656  300.10   <2e-16
## hour_groupsevening:rateCar_LaggedHour      0.572480   0.003810  150.26   <2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.881992   0.002559  344.71   <2e-16
## hour_groupsnight:rateCar_LaggedHour        0.563815   0.011286   49.96   <2e-16
## hour_groupsnoon:rateCar_LaggedHour         0.928103   0.002810  330.23   <2e-16
##                                               
## (Intercept)                                ***
## hour_groups8heures:rateCar_LaggedHour      ***
## hour_groups9heures:rateCar_LaggedHour      ***
## hour_groupsafternoon:rateCar_LaggedHour    ***
## hour_groupsevening:rateCar_LaggedHour      ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour        ***
## hour_groupsnoon:rateCar_LaggedHour         ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.421 on 35052 degrees of freedom
## Multiple R-squared:  0.8753, Adjusted R-squared:  0.8753 
## F-statistic: 3.516e+04 on 7 and 35052 DF,  p-value: < 2.2e-16

Toutes les variables ont une contribution significative d’après les tests de student. Calculons le score de validation croisée pour ce modèle simple pour obtenir

CV_gam = function(formula, df_train, no_folds){
  K = no_folds 
  data = df_train
  folds <- cut(seq(1,nrow(data)),breaks=K,labels=FALSE)
  fold_scores = c()
  start_time = Sys.time()
  for(i in 1:K){
    testIndexes <- which(folds==i,arr.ind=TRUE)
    testData <- data[testIndexes, ]
    trainData <- data[-testIndexes, ]
    g = gam(data=trainData, formula=formula)
    Y_test = testData$rateCar
    Y_predict = predict(g, newdata=subset(testData, select = -c(rateCar)))
    fold_scores = c(fold_scores, rmse(Y_test, Y_predict))
  }
  exec_time = Sys.time()-start_time
  score = mean(fold_scores)
  ret = data.frame(score = score, no_folds = no_folds, exec_time = exec_time, 
                   formula = paste(as.character(formula)[c(2,1,3)], collapse = " "))
  return(ret)
}
knitr::kable(CV_gam(form1, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1)Model score from CV

score	no_folds	exec_time	formula
2.42125	8	1.119009 secs	rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour

Afin de choisir la prochaine variable à rajouter pour expliquer rateCar, nous les explorons une par une :

1. hour : Circulaire, essayons différents k.

form1_hour = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(hour, bs="cc", k = 24, by=weekendsIndicator)
mod1_hour = gam(form1_hour, data=df_train)
summary(mod1_hour)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + 
##     s(hour, bs = "cc", k = 24, by = weekendsIndicator)
## 
## Parametric coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                1.101257   0.027548   39.98  < 2e-16
## hour_groups7heures:rateCar_LaggedHour      0.711020   0.135167    5.26 1.45e-07
## hour_groups8heures:rateCar_LaggedHour      2.304156   0.107999   21.34  < 2e-16
## hour_groups9heures:rateCar_LaggedHour      2.343870   0.033249   70.49  < 2e-16
## hour_groupsafternoon:rateCar_LaggedHour    1.026931   0.007463  137.61  < 2e-16
## hour_groupsevening:rateCar_LaggedHour      0.525550   0.006976   75.33  < 2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.819606   0.004874  168.18  < 2e-16
## hour_groupsnight:rateCar_LaggedHour        0.655250   0.016846   38.90  < 2e-16
## hour_groupsnoon:rateCar_LaggedHour         0.867199   0.004684  185.13  < 2e-16
##                                               
## (Intercept)                                ***
## hour_groups7heures:rateCar_LaggedHour      ***
## hour_groups8heures:rateCar_LaggedHour      ***
## hour_groups9heures:rateCar_LaggedHour      ***
## hour_groupsafternoon:rateCar_LaggedHour    ***
## hour_groupsevening:rateCar_LaggedHour      ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour        ***
## hour_groupsnoon:rateCar_LaggedHour         ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                              edf Ref.df      F p-value    
## s(hour):weekendsIndicator0 21.94     22 203.63  <2e-16 ***
## s(hour):weekendsIndicator1 21.44     22  75.37  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.894   Deviance explained = 89.4%
## GCV = 5.0041  Scale est. = 4.9966    n = 35060

knitr::kable(CV_gam(form1_hour, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_hour)Model score from CV

score	no_folds	exec_time	formula
2.2425	8	5.049379 secs	rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(hour, bs = “cc,” k = 24, by = weekendsIndicator)

#test_score = rmse(predict(mod1_hour, newdata = df_test),  df_test$rateCar)
#print(test_score)

2. weekdays :

form1_weekdays = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(weekdays, bs="cc", k = 4,                                  by=as.factor(winterHolidaysIndicator)) 
# + s(weekdays, bs="cc", by=summerHolidaysIndicator) + s(weekdays, bs="cc", by=bankHolidaysIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
mod1_weekdays = gam(form1_weekdays, data=df_train)
summary(mod1_weekdays)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + 
##     s(weekdays, bs = "cc", k = 4, by = as.factor(winterHolidaysIndicator))
## 
## Parametric coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                0.859531   0.025096   34.25   <2e-16
## hour_groups7heures:rateCar_LaggedHour      0.853540   0.074459   11.46   <2e-16
## hour_groups8heures:rateCar_LaggedHour      2.186405   0.039039   56.01   <2e-16
## hour_groups9heures:rateCar_LaggedHour      1.928966   0.013898  138.80   <2e-16
## hour_groupsafternoon:rateCar_LaggedHour    1.095019   0.003699  296.01   <2e-16
## hour_groupsevening:rateCar_LaggedHour      0.569527   0.003870  147.16   <2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.880178   0.002603  338.14   <2e-16
## hour_groupsnight:rateCar_LaggedHour        0.555410   0.011345   48.96   <2e-16
## hour_groupsnoon:rateCar_LaggedHour         0.925539   0.002894  319.86   <2e-16
##                                               
## (Intercept)                                ***
## hour_groups7heures:rateCar_LaggedHour      ***
## hour_groups8heures:rateCar_LaggedHour      ***
## hour_groups9heures:rateCar_LaggedHour      ***
## hour_groupsafternoon:rateCar_LaggedHour    ***
## hour_groupsevening:rateCar_LaggedHour      ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour        ***
## hour_groupsnoon:rateCar_LaggedHour         ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                                       edf Ref.df     F p-value
## s(weekdays):as.factor(winterHolidaysIndicator)0 1.948e+00      2 46.78  <2e-16
## s(weekdays):as.factor(winterHolidaysIndicator)1 5.889e-09      2  0.00   0.532
##                                                    
## s(weekdays):as.factor(winterHolidaysIndicator)0 ***
## s(weekdays):as.factor(winterHolidaysIndicator)1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.876   Deviance explained = 87.6%
## GCV = 5.8457  Scale est. = 5.8439    n = 35060

#plot(mod1_weekdays, residuals=T, rug=T, se=F, pch=20)
knitr::kable(CV_gam(form1_weekdays, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_wd)Model score from CV

score	no_folds	exec_time	formula
2.415	8	1.290045 secs	rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(weekdays, bs = “cc,” k = 4, by = as.factor(winterHolidaysIndicator))

3. state :

form1_state = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(state, bs="cr", k = 5)
# + s(weekdays, bs="cc", by=summerHolidaysIndicator) + s(weekdays, bs="cc", by=bankHolidaysIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
mod1_state = gam(form1_state, data=df_train)
#summary(mod1_state)
#plot(mod1_weekdays, residuals=T, rug=T, se=F, pch=20)
knitr::kable(CV_gam(form1_state, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_state)Model score from CV

score	no_folds	exec_time	formula
2.42125	8	0.750025 secs	rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(state, bs = “cr,” k = 5)

3. toy (time of year):

form1_toy = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(toy, bs="cr", k = 6)
# + s(weekdays, bs="cc", by=summerHolidaysIndicator) + s(weekdays, bs="cc", by=bankHolidaysIndicator)
# mod1_hour = gam(form1_hour, data=df_train)
# summary(mod1_hour)
mod1_toy = gam(form1_toy, data=df_train)
summary(mod1_toy)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + 
##     s(toy, bs = "cr", k = 6)
## 
## Parametric coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                0.831406   0.024825   33.49   <2e-16
## hour_groups7heures:rateCar_LaggedHour      0.882487   0.074472   11.85   <2e-16
## hour_groups8heures:rateCar_LaggedHour      2.209279   0.038941   56.73   <2e-16
## hour_groups9heures:rateCar_LaggedHour      1.937016   0.013848  139.88   <2e-16
## hour_groupsafternoon:rateCar_LaggedHour    1.097766   0.003682  298.14   <2e-16
## hour_groupsevening:rateCar_LaggedHour      0.573190   0.003842  149.18   <2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.882456   0.002584  341.51   <2e-16
## hour_groupsnight:rateCar_LaggedHour        0.563128   0.011324   49.73   <2e-16
## hour_groupsnoon:rateCar_LaggedHour         0.928841   0.002852  325.67   <2e-16
##                                               
## (Intercept)                                ***
## hour_groups7heures:rateCar_LaggedHour      ***
## hour_groups8heures:rateCar_LaggedHour      ***
## hour_groups9heures:rateCar_LaggedHour      ***
## hour_groupsafternoon:rateCar_LaggedHour    ***
## hour_groupsevening:rateCar_LaggedHour      ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour        ***
## hour_groupsnoon:rateCar_LaggedHour         ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value    
## s(toy) 4.676  4.948 4.658 0.000219 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.875   Deviance explained = 87.5%
## GCV = 5.8582  Scale est. = 5.8559    n = 35060

#plot(mod1_weekdays, residuals=T, rug=T, se=F, pch=20)
knitr::kable(CV_gam(form1_toy, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_toy)Model score from CV

score	no_folds	exec_time	formula
2.42125	8	0.7370381 secs	rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(toy, bs = “cr,” k = 6)

3. nbCar_LaggedHour:

for (h in 0:23){
  g = ggplot(data=df_train[(df_train$hour == h),], aes(y = rateCar, x = nbCar_LaggedHour))+ 
  geom_point() + xlab("nbCar_laggedHour") +
  ylab("rateCar (présent)") + 
  labs(title = paste0("Heure ", h))
plot(g)
}

form1_nbCarHour = rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(nbCar_LaggedHour, bs="cr", k = 6, by=hour_groups)
mod1_nbCarHour = gam(form1_nbCarHour, data=df_train)
summary(mod1_nbCarHour)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + 
##     s(nbCar_LaggedHour, bs = "cr", k = 6, by = hour_groups)
## 
## Parametric coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                1.845837   0.082954  22.251  < 2e-16
## hour_groups7heures:rateCar_LaggedHour      0.905683   0.202185   4.479 7.51e-06
## hour_groups8heures:rateCar_LaggedHour      1.560141   0.162439   9.604  < 2e-16
## hour_groups9heures:rateCar_LaggedHour      1.361042   0.066729  20.396  < 2e-16
## hour_groupsafternoon:rateCar_LaggedHour    0.929243   0.009734  95.468  < 2e-16
## hour_groupsevening:rateCar_LaggedHour      0.536947   0.010040  53.481  < 2e-16
## hour_groupsevening_rush:rateCar_LaggedHour 0.818969   0.005512 148.568  < 2e-16
## hour_groupsnight:rateCar_LaggedHour        0.455403   0.025760  17.679  < 2e-16
## hour_groupsnoon:rateCar_LaggedHour         0.804931   0.006956 115.710  < 2e-16
##                                               
## (Intercept)                                ***
## hour_groups7heures:rateCar_LaggedHour      ***
## hour_groups8heures:rateCar_LaggedHour      ***
## hour_groups9heures:rateCar_LaggedHour      ***
## hour_groupsafternoon:rateCar_LaggedHour    ***
## hour_groupsevening:rateCar_LaggedHour      ***
## hour_groupsevening_rush:rateCar_LaggedHour ***
## hour_groupsnight:rateCar_LaggedHour        ***
## hour_groupsnoon:rateCar_LaggedHour         ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                               edf Ref.df      F p-value    
## s(nbCar_LaggedHour):hour_groups7heures      3.610  3.884  10.44  <2e-16 ***
## s(nbCar_LaggedHour):hour_groups8heures      4.083  4.374  95.04  <2e-16 ***
## s(nbCar_LaggedHour):hour_groups9heures      4.086  4.277 159.35  <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsafternoon    4.985  5.000  36.38  <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsevening      4.847  4.981  54.04  <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsevening_rush 4.761  4.949  17.89  <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsnight        3.250  3.682 132.29  <2e-16 ***
## s(nbCar_LaggedHour):hour_groupsnoon         4.797  4.979  27.27  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.885   Deviance explained = 88.5%
## GCV = 5.4364  Scale est. = 5.4297    n = 35060

#plot(mod1_nbCarHour, residuals=T, rug=T, se=T, pch=20)
#draw(mod1_nbCarHour, residuals=T)
knitr::kable(CV_gam(form1_nbCarHour, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_nbCarHour)Model score from CV

score	no_folds	exec_time	formula
2.4525	8	3.326821 secs	rateCar ~ hour_groups * rateCar_LaggedHour - hour_groups - rateCar_LaggedHour + s(nbCar_LaggedHour, bs = “cr,” k = 6, by = hour_groups)

3.3 Arbres aléatoires

Après cet exercise de bricolage avec GAM, nous passons à d’autres modèles dont le tuning pourra se faire de manière plus automatique. Pour commencer, nous nous intéressons à des arbres CART et nous nous servons du package rpart pour les implémenter. Au lieu d’examiner la contribution de chaque variable, le choix des \(k\) etc., il suffit de tuner quelques paramètres de l’algorithme. Vu que nous ferons un tuning étendu pour les forêts aléatoires après, nous nous focalisons ici sur un seul paramètre qui est de plus haute importance : le “complexity parameter” cp. Dans la construction de l’arbre, si le gain de performance après une découpe n’est pas meilleur d’un facteur de cp, alors la découpe n’est pas réalisée. Plus cp est grand, moins l’arbre sera alors complexe, dont le nom.

Pour optimiser cp, nous utilisons de la validation croisée et une recherche sur deux grilles : la première de 0.0 à 0.1 et la seconde, plus fine, de 0.0 à 0.01. Nous affichons les scores de CV obtenues (voir le fichier 1_heure/tree.R pour le code) :

Figure 3.3: Recherche du paramètre cp optimal

La valeur optimale d’après cette CV est alors cp=0, c’est-à-dire que l’on accepte toutes les découpes et que notre arbre est un arbre profond. On maintient à leur valeur par défaut les autres paramètres de CART et on entraîne le modèle sur notre train dataset. Nous faisons ceci dans deux cas différents qui représentent deux dimensions locales différentes : rappelons-nous que pour la complétion, nous avions rajouté des informations sur les voisins pour enrichir les données. Pour voir si ces mêmes informations nous servent aussi maintenant, nous entraînons un arbre avec et un CART sans données voisines. Voici les scores sur le test set dans les deux cas :

Table 3.2: Performances des arbres aléatoires

Arbre avec voisins	Arbre sans voisins
3.783	3.582

Premièrement, on note dans les deux cas une performance légèrement meilleure que le modèle naïf qui recopie l’heure d’avant. Ce qui est pourtant frappant est que le rajout d’information sur les voisins baisse le score, indiquant qu’il y a soit de l’overfitting, soit un problème de grande dimension qui perturbe le bon fonctionnement des arbres aléatoires. En tout cas, ce résultat sème le doute sur l’utilité de la composante spatiale de nos données. Mais peut-être, il faut juste utiliser des techniques plus puissantes, comme les forêts aléatoires qui rassemblent plusieurs arbres en utilisant du bagging.

3.3.1 Forêts aléatoires

Pour cette méthode d’ensemble basé sur les arbres de Breimann, le package ranger nous servira à construire les forêts et le package vip aidera à visualiser l’importance des variables. Comme pour les arbres, il y a des paramètres à optimiser avec de la CV pour que la forêt apprenne bien mais pas trop la circulation parisienne. Nous nous focalisons sur deux paramètres que nous optimisons l’un après l’autre. C’est une approche moins prometteuse qu’une recherche sur grille, mais elle diminue le temps de calcul et nous permet de visualiser la recherche de paramètres :

• Le premier paramètre à optimiser est min.node.size (minimal node size). L’algorithme ne découpe pas une feuille dont le nombre d’éléments est inférieur ou égal à sa valeur. Comme cp avant, ce paramètre contrôle la taille des arbres.

Figure 3.4: Recherche du paramètre min.node.size optimal

On remarque que l’erreur décroît jusque la valeur 1 avec un pic négatif à 9. Choisir min.node.size=1 donnerait des arbres de profondeur maximale, rendant le modèle susceptible à sur-apprendre. Par conséquent nous nous contentons de choisir 9, malgré les fluctuations importantes des CV-scores dans le graphique à droite. Ces flucuations pourraient être dues au fait qu’ici on a fait une validation croisée habituelle (voir nos explications sur CV en haut) qui mélange les lignes corrélées avant de les séparer en blocs. Dans un cadre iid, nous nous attendrions à la forme typique d’un “coude” par exemple.

• Le deuxième paramètre est mtry : le nombre de variables, tirées au hasard, qui sont considérées pour la construction de chaque arbre. Rendre ce paramètre plus petit réduit alors la corrélation entre les arbres de la forêt. Pour l’optimiser, nous utilisons deux de nos approches CV :

Figure 3.5: Recherche du paramètre mtry optimal

Celle-ci est l’unique fois où nous affichons la méthode de validation croisée par timeslice. On remarque que l’allure globale des courbes sont identiques mais que le coude est plus lisse avec timeslice. On choisit de poser \(\textbf{mtry} = 10\), qui semble être un bon choix d’après les deux graphiques.

Le dernier paramètre à déterminer est ntree, le nombre d’arbres dans la forêts. Ce n’est pas un paramètre à optimiser comme les 2 précédents car le sur-apprentissage dans les forêts aléatoires n’est pas lié au nombre d’arbres. On optimise ntree afin de déterminer à partir de quel nombre d’arbre on gaspille du temps d’exécution pour un faible gain de performance.

Figure 3.6: Recherche du paramètre ntree ‘’optimal’’

On voit que le temps d’exécution est linéaire selon le nombre d’arbres et que la courbe de l’erreur marque un coude à partir de \(150\), valeur que l’on choisit pour **ntree*.

Avec ces trois paramètres choisis, nous entraînons la forêt aléatoire et l’évaluons sur le test set. Comme avec les arbres, nous faisons ceci deux fois, une fois avec voisins et une fois sans.

Table 3.3: Performances des forêts aléatoires

Forêt aléatoire avec voisins	Forêt aléatoire sans voisins
3.092	3.129

Comme attendu, les performances des forêts aléatoires sont généralement meilleures que celles des arbres aléatoires. Pourtant, contrairement à ces derniers, les arêtes voisins semblent avoir un effet plutôt positif. Si nous nous intéressons davantage au rôle des voisins, nous pouvons exploiter le fait que les forêts aléatoires livrent des scores d’importance pour chaque variable. Sur l’arête pont amont - pont austerlitz, par exemple, le ranking d’importance des variables est le suivant :

Figure 3.7: Scores d’importance dans l’ordre décroissante (15 premières arêtes) pour l’arête ‘’pont amont - pont austerlitz’’

Sur cette arête, les voisins dominent même les scores d’importance. Ceci nous donne une motivation pour davantage explorer le rôle des voisins et ainsi le rôle de la dimension spatiale du problème de prédiction de la circulation :

Une question immédiate est quelles arêtes profitent le plus des données de leurs voisins. Nous essayons de visualiser ça sur le graphique suivant où pour chaque arête du graphe routier, nous visualisons la différence des erreurs test avec et sans voisins. Le trait vertical noir pointillé marque la séparation entre les rues à l’intérieur Paris et les parties du périphérique. Ces dernières se trouvent à droite de la ligne. Le trait horizontal rouge pointillé est un seuil arbitraire à partir duquel nous interprétons que l’ajout des voisins impacte la performance de la forêt significativement. Les arêtes qui dépassent le seuil en bas (et souffrent donc de l’inclusion de leurs voisins) sont représentées par une croix rouge, les autres par une croix noire.

TODO Jakob discuter avec Guillaume. Seules 12 forêts sur 69 ont une performance ‘’significativement’’ meilleure lorsqu’on ajoute les voisins. De plus, pour les arêtes du périphérique, on remarque que les différences sont plus importantes que pour Paris, qu’inclure les voisins est néfastes (à 2 exceptions près) et qu’en général, la circulation sur ces arêtes est moins bien prédite. En effet, les scores pour Paris sont 2.747 (avec voisins) et 2.898 (sans voisins), et les scores du périphérique sont 4.069 (avec voisins) et 3.787 (sans voisins).

Conclusion

“Covid Index.” visité le 26.12.2021. https://www.bsg.ox.ac.uk/research/research-projects/covid-19-government-response-tracker.

“Données Du Comptage Routier.” visité le 26.12.2021. https://parisdata.opendatasoft.com/explore/dataset/comptages-routiers-permanents-historique/information/.

“Jours Fériés.” visité le 26.12.2021. https://www.data.gouv.fr/en/datasets/jours-feries-en-france/.

“Les Archives Du Calendrier Scolaire.” visité le 26.12.2021. https://www.education.gouv.fr/les-archives-du-calendrier-scolaire-12449.

“Observations Météorologiques Historiques En France.” visité le 26.12.2021. https://public.opendatasoft.com/explore/dataset/donnees-synop-essentielles-omm/table/?sort=date&q=%C3%AEle+de+france&refine.nom=ORLY&q.timerange.date=date:%5B2013-12-31T23:00:00Z+TO+2021-01-01T22:59:59Z%5D.

“Paris Data.” visité le 26.12.2021. https://parisdata.opendatasoft.com/page/home/.

“Progressive Cross-Validation.” visité le 22.02.2022. https://medium.com/@soumyachess1496/cross-validation-in-time-series-566ae4981ce4.

Van Buuren, Stef. 2018. Flexible Imputation of Missing Data. CRC press.